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Abstract. There are two types of two-photon transitions in molecular aggregates, 
that is, non-local excitations of two monomers and local double excitations to some 
higher excited intra-monomcr electronic state. As a consequence of the inter-monomer 
Coulomb interaction these different excitation states are coupled to each other. 
Higher excited intra-monomer states are rather short-lived due to efficient internal 
conversion of electronic into vibrational energy. Combining both processes leads to 
the annihilation of an electronic excitation state, which is a major loss channel for 
establishing high excitation densities in molecular aggregates. 

Applying theoretical pulse optimization techniques to a Frenkel exciton model it is 
shown that the dynamics of two-exciton states in linear aggregates (dimer to tetramer) 
can be influenced by ultrafast shaped laser pulses. In particular, it is studied to 
what extent the decay of the two-exciton population by inter-band transitions can be 
transiently suppressed. Intra-band dynamics is described by a dissipative hierarchy 
equation approach, which takes into account strong exciton-vibrational coupling in 
the non-Mar kovian regime. 
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1. Introduction 

Molecular aggregates continue to provide inspiration and challenges to experiment and 
theory [JJ El E]- In terms of a microscopic understanding of the dynamics of elementary 
excitations recent advances due to time-resolved nonlinear spectroscopy have been 
tremendous [U [5] . However, despite the success of laser pulse control in rather different 
areas of research [6], relatively little emphasis has been put on exciton dynamics. On the 
experimental side, feedback laser control has been applied to manipulate the branching 
ratio between internal conversion (IC) and energy flow in the light-harvesting antenna of 
purple bacteria [7J . There has been a number of simulations for light-harvesting systems 
by May and Bruggemann et al. where the focus was put on the transient generation of 
a single localized excitation within the aggregate 01211101111]; non-biological aggregates 
have received no attention so far. 

The quantum dynamics of excitons is well founded within the Frenkel exciton 
approach, which starts from a classification of aggregate's excitation states in terms 
of the number of simultaneously excited monomers [T2l IT3] . Under weak irradiation 
conditions only a single excitation will be present, but in nonlinear optical experiments 
or in the presence of strong irradiation multiple excited states play an important role. 
The organic building blocks of molecular aggregates have a multitude of electronically 
excited states, i.e. there is always a state S n such that the So-Si excitation energy 
approximately matches that of a Si-S n transition. As a consequence one needs to 
distinguish between local double excitations (LDE) and nonlocal double excitations 
(NDE) where the two excitations are localized at different monomers. Note that the 
NDE should not be confused with bi-exciton states, i.e. bound states formed by two 
excitons in the presence of different permanent dipoles in the involved electronic states 
[13] . LDE and NDE can couple via the Coulomb interaction. The manifestation of 
this coupling in nonlinear spectroscopy has been investigated in Refs. [HI [151 EH]- The 
presence of LDEs has important consequences for the exciton dynamics since it leads to 
exciton-exciton annihilation (EEA) by virtue of an intramolecular IC triggered by non- 
adiabatic electronic transitions. The presence of this process in molecular aggregates is 
well established [171 [TBI US] and various theoretical approaches exist [201 ED E21 E31 El] , 
although it must be stated that a first principles determination of the respective IC 
rates is still out of reach. 

Understanding exciton dynamics in aggregates is impossible without taking into 
account the effect of exciton- vibrational coupling [31 [121 ES]- Here one can distinguish 
between approaches which either account for all vibrations on the same footing, i.e. in 
terms of a heat bath, or treat a few active vibrational degrees of freedom explicitly, 
but coupled to the heat bath of the remaining modes. Needless to say, the latter 
approach is more demanding as the dimension of the density matrix increases rapidly 
with the number of explicit modes. Early investigations along these lines therefore have 
been restricted to molecular dimer models with one active vibrational coordinate per 
monomer [261 EJJ EH E2] ■ More recent approaches include a Green's function description 
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of exciton dynamics [30], a mult i- configuration time- dependent Hartree simulation [31] 
or the Frenkel excitonic polaron treatment [32J. These methods cannot treat finite 
temperature effects due to the coupling to a further heat bath in a microscopic manner. 

The coupling of the electronic degrees of freedom to a single or site-specific heat bath 
is usually treated using dissipation theory [3] . Recently, the so-called hierarchy equation 
(HE) method [33], EH [35j EEl EH [381 EH] has enjoyed great popularity as it promises a non- 
perturbative and non-Markovian description of exciton dynamics, which is numerically 
equivalent, e.g. to the influence functional approach within a path integral formulation 
[3D]. In the present context one should note that a distinction between different 
types of modes, i.e. low-frequency solvent modes versus high-frequency intramolecular 
vibrations, can be introduced via the spectral density. This is of particular relevance 
since the chromophores used to assemble artificial molecular aggregates usually show 
vibronic side bands in their absorption spectra pointing to the prominent role of 
intramolecular high-frequency vibrations |41j. 

The focus of the present contribution is on the laser control of dissipative two- 
exciton dynamics in linear aggregates. Thereby we take into account the effect of EEA 
and of a coupling to a heat bath being composed of an effective high-frequency mode 
as well as of a continuous distribution of low-frequency modes. For the solution of the 
density matrix equations of motion we will use a HE approach. Two-exciton populations 
are usually quenched by EEA. Since the latter is a local process we will ask the question 
whether an excitonic wave packet can be prepared such as to transiently suppress EEA 
by its composition in terms of LDE states. In the following Section [2] we will outline 
the density matrix approach and its interface to optimal control theory (OCT). Section 
[3] will start with a discussion of the field free dynamics focussing on the comparison 
between the Markovian and non-Markovian limits. Subsequently, laser driven dynamics 
is discussed. The results are summarized in Sec. HI 



2. Theory 

2.1. Frenkel Exciton Hamiltonian 

Consider an aggregate composed of N monomers, each carrying three adiabatic 
electronic states (a = g,e,f) corresponding to the So, Si, and some S n state. Thus we 
have the adiabatic electronic basis of local states \m, a), m = 1, . . . , N. The electronic 
states of the aggregate can be classified as the zero excitation (ground) state [3] 

|0> = n m |m,s>, (1) 
the singly-excited state 

\m) = \m,e) Ii n ^ m \n,g) , (2) 

the doubly-excited (LDE and NDE) states 

\mm) = \m,f)U n ^ m \n,g), (3) 
\mn) = \m,e) \n,e) n fe ^ m>n \k,g) . (4) 
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Restricting ourselves to these excitation the completeness relation reads 

1 = |0) (0| + y~^(|m) (m| + \mm) (mm\) + - \mn) (mn\ . (5) 

The electronic Hamiltonian can be written as 

H c \ = H cx + if fic i d (t) + Hie , (6) 
with the bare exciton Hamiltonian, 

£ex = flj? + fl£> + flg> , (7) 

being composed of the ground state part 

H2? = Eo |0) (0| , (8) 
where E is the electronic ground state energy, the single-exciton Hamiltonian 

Hg> = E m,e \m) (m\ + J mn H (n\ , (9) 

m,e mj^n 

and the two-exciton Hamiltonian 

= Y Em 'f l mm ^ ^ mm l + o X^ m > e + ^"-^ l mn ^ ^ mn l 

+ J n fc |mn) (mA;| + (j^ |mn) (mm| + h.c.) . (10) 

Here, E m>e and E m j are the energies of electronic excitation at site m for So-Si and 
S -S n , respectively. Further, J mn is the coupling between site m and site n leading to 
single exciton transfer and Jmn is the Coulomb coupling responsible for the two-exciton 
dynamics. 

It is customary to introduce Frenkel exciton eigenstates which follow from the 
diagonalization of the exciton Hamiltonian 

H cx \a) = E a \a) . (11) 

The eigenstates can be decomposed in terms of local excitation states, \a) (a = 
0, m, mm, mn), as follows 

l«> = $^C a («)|a). (12) 

a 

The eigenstates separate into .M-exciton manifolds (here M. = 0,1,2). Frequently, 
we will also use a notation where the states in the .M-exciton manifolds are counted 
according to increasing energy, i.e. \A4k)- 

In Eq. (g) the aggregate is assumed to interact with an external laser field treated 
in dipole approximation, i.e. 

H Aeld (t) = -E(t)jl, (13) 

where the field E(t) is oriented parallel to the transition dipoles. In the local Frenkel 
exciton basis the dipole operator for the aggregate reads 



A = ^ I AW \m) (0| + fi mJ \mm) (ra| + ^ AV 

m \ n=£m 



\mn) (m\ 1 +h.c. (14) 
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Here, n m ,e{y>m,f) is the transition dipole for So-Si(Si-S n ) transition of monomer m. 

Finally, we consider the IC process described by Hie in Eq. (JHfc. IC has its 
origin in the breakdown of the Born-Oppenheimer approximation, e.g. at conical 
intersections. The respective non-adiabaticity operator triggering transitions between 
adiabatic electronic states is proportional to the momentum operator P m ^ of the involved 
nuclear degrees of freedom (coordinates Q m ^) counted by the index £ at monomer m. 
The non-radiative life time of the Si state is usually in the nanosecond range. Therefore 
we will restrict ourselves to the IC between adiabatic states \m, /) and \m,e). Hence 
the coupling becomes 

Hie ^E^Efi^, 

n^ -* = \m) (mm\ + \mm) (m\ . (15) 
In the following we will assume that the coupling strength is site-independent, i.e. 

JIC) (IC) 

2.2. Exciton- Vibrational Coupling: System-Bath Model 

Exciton-vibrational coupling leading to intra-band phase and energy relaxation is 
introduced in the spirit of the system-bath approach, i.e. we have 

H = Hs + Hb + Hs-b 

= H c \ + H vih + H c \^ vih . (16) 

The bath modes are treated in harmonic approximation, i.e. Hvih is the standard 
harmonic oscillator Hamiltonian. Concerning the system-bath coupling we will 
distinguish between two types of modes (see also Ref. [12] )• First, local high-frequency 
modes, {q m ,£,}, which usually give rise to a vibrational progression in the monomer 
absorption spectrum [UJ. Second, global low-frequency modes, {x^}, contributing via 
a continuous spectrum. Hence we have 

^ci-vib= ^^ H) + i^ (L) . (17) 

m 

For the coupling to the local high-frequency modes we will use the model Hamiltonian 

^m H) = $^<W 9m\, e \ m ) H + gSh \ mm ) ( mm \ 

+ 9 ^l,e \ mn ) ( mU \ ■ ( 18 ) 

n^ra 

In the following we will assume that the coupling is the same for all monomers, i.e. 
9m\a = dfa- Further, we will use = g^} and g^j = ng^ [13]. With these 



assumptions Eq. (18) can be written as 

h m — \m) (m\ + K \mm) (mm\ + \mn) (mn\ . (19) 
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For the coupling to the low-frequency modes we will invoke the same approximations 
and arrive at the Hamiltonian 

#W = /,(tot) ^(L) % 

fe (tot) = ^L. (20) 

m 

2.3. Dissipation Models 

2.3.1. Internal Conversion The IC rate between S„ and Si states is rather large for 
typical chromophores, reaching time scales on the order of about 100 fs jB]. On this 
time scale the S n population decays into the Si state, thereby passing many electronic 
states (for the case of perylene bisimides n would be of the order of about 30 [2]). 
Hence it is justified to assume that the memory time associated with the IC process is 
even shorter than 100 fs and EEA can be treated in Markovian approximation. Treating 



Hie in Eq. ( 15 ) as a perturbation of the bare exciton Hamiltonian within second order 
perturbation theory one can write the contribution to the Quantum Master Equation 
for the reduced exciton density operator p in terms of the Lie operator lZ <ylc " > which 
operates on p as [3] 

ifSRP^p = Yl [n£ C) , (S m p - 0j] ■ (21) 



Here E m is defined as 



POO 

/ dta^ c \t)exp(-iH e J/h)fl^ c) exp(iH ex t/h) (22) 
Jo 



and «£ C '(t) is correlation function of the part ^2^g^P m ^ °f Eq. (15) (note that we 
assume that the momenta at different sites are uncorrelated) . By using the completeness 
relation for the eigenstates, one has 

s m = E nJS» J £ c) ((£/» - E «)/ n ) l«> W • C 23 ) 

Here, J^\oj) = dtexp(iojt)am (t) is the Fourier transformation of the correlation 

(IC) 

m,a/3 



function and n£ c 2 5 = (a| \/3). For simplicity the Debye spectra is used below 



jg C \u)= mc u J -^- r (24) 

2.3.2. Separation of Time Scales in the Response Function of the Oscillator Bath In 
general the influence of the bath degrees of freedom due to exciton-vibrational coupling 
cannot be treated in Markovian approximation. It is fully characterized by the response 
function, a(t), of the bath, which in turn is determined by the spectral density function 
J(u) 

1 [°° 

a(t) =— did J(lo) [coth(/3faj/2) cos(a;£) — i sm(cut)\ , 
^ Jo 
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J ^) = lJ2— (25) 

Here, /i^, and are the coupling constant, reduced mass, and frequency, respectively, 
of the £th oscillator and /3 = l/k^T. In the following we will specify the spectral density 



to the coupling models discussed in Section 2.2 



First, we consider the case of the local high-frequency mode, Eq. (18) for which 
we assume a damped Brownian oscillator model to hold. This can be described by the 
following spectral density 



J hM = -i— 2 , 2 . (26) 

Below we will assume the underdamped limit, where the central mode frequency Uq is 
much larger than the cutoff Ah- For this spectral density the response function can be 
expanded into a sum of exponentials as follows 

a H (t) = {-exp(^) [coth(^/30 x /2) - 1] + exp(fi 2 t) [coth(^/3fi 2 /2) - 1]} 

-2 W 3 A^/?£ 



, y k + uif - a^| 



AT; 



II 



= he~ nkt + c i H) exp(-^t) + a<t mt) (t). (27) 

k=l k=l 

where fl 1)2 = A H /2 ± if, f = ^/wg _ A|/i. 



In Eq. (27), z/^ = 2irk/(h/3)(k > 0) is the fcth Matsubara frequency of the bath, and 
Ah is the smallest integer satisfying z^v H +i 3> Ah- As discussed in Ref. [33] the response 
function can be split into two parts: a long memory contribution (the first two terms in 
the last line of Eq. (27)) and a short memory part OTt \t). The latter one contains 
all the terms with short memory in the Matsubara summation and will be treated in 
Markovian approximation within the hierarchy equations to be discussed below. 

The low-frequency bath will be described by the Debye spectral density 



j 2 + A 2 l 

where t]l is the coupling strength and Al is the frequency cutoff of the bath. The 



response function in Eq. (25) then can be expressed as follows 



a L (t) = [cot imL /2) - i] + £ |^ 



2 

k=l 

N L 

E 

k=0 



ci L) exp(-^) + 4 short) (t), (29) 



where u = A L , v k = l-akj {h(3){k > 0) is the fcth Matsubara frequency of the bath, and 
Al is the smallest integer satisfying vn^+i S> Al. The same time scale separation has 
been introduced as for the high-frequency bath. 
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2.4- Hierarchy Equations for the Dissipative Multi-exciton Dynamics 



In order to describe the non-Markovian and non-perturbative exciton dynamics we will 
utilize a HE approach to propagate the reduced exciton density matrix in eigenstate 
representation. Specifically we have employed the stochastic decoupling procedure due 
to Shao and coworkers which is briefly sketched in the Appendix [37J EH] ■ Since the HE 
approach starts from the influence functional [33], 06], it is based on the fact that the 
response function for bath is written as a sum of exponentials (see previous section). 

The resulting HEs of motion for the dissipative exciton dynamics of the aggregate 
are given by 

N 2 N Nn N L 



d 



vm=l j=l 



m=l j=l 



k=0 



+ 



H cx + H Rcld (t), PABV (t)\ +J2^ 1 \/(V j + l)\cf ) \ U tot \p ABV+ (t) 



3=0 




(t) 



m >PAB+ 3 V 



(t) 



(t) 



— 01 « 



h 



(t) 



(t) 




h {t0t) ,p ABVr (t) 



ihcf}{h^\p ABVr {t) 



+ c (L) 



h (tot) , P ABV -(t) 



(30) 



where A is an iV x 2 integer matrix, with its matrix element A m j corresponding the 
contribution from Qj in m-th molecule (Eq. (27)). B is an x A^h integer matrix with 
its matrix element B m j counting the j-th Matsubara frequency in the m-th molecule 
(Eq. (27)). And V is an ]Vl + 1 dimensional integer vector. Vo refers to the cutoff 
frequency and Vj (j ^ 0) to the j-th Matsubara frequency in the low frequency bath 
(Eq. (29)). The operator " + / — " is defined as [M mJ -] = M n i ± 5 mn 5ji for the matrix 
and [Vj ] = 14 ± 5kj for the vector. 

VJ is the Redfield super-operator accounting for the exciton-exciton annihilation 
and the short memory effects of the high and low frequency environments, [•, •] denotes 
the commutator and {•, •} denotes the anti-commutator. For the definition of otk,± an d 
Appendix. 

In the above equations, the first term PQOO ls ^ e re duced density matrix for the 
exciton and the other elements reflect the finite memory effect due to the exciton- 
vibration coupling. The initial condition for the hierarchy is Pqqq(0) = p(O) and 
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f>ABv{fy — if any of the elements of A, B and V is nonzero. 

In the Markovian limit, the HEs are reduced to the Quantum Master Equation 



lh dt P 



P 



E(t) \fi, p] — ihlZp, 



(31) 



where 71 is the Redfield superoperator 



N 



(32) 



Here TZ^ is the superoperator accounting for the IC as defined in Eq. (21). H( tot ) 



and S m are the dissipative operators defined in the similar way as Eq. (22) for the low 
frequency and high frequency bath in m-th molecule, respectively. 



2.5. Optimal Control Theory 

For the design of laser fields for driving the exciton dynamics we will employ optimal 
control theory (see, eq. Ref. |47J; the present implementation follows Ref. [3]). Here 
the goal is to find a laser field E(t) such as to optimize a target at a certain time t = T. 
This can be cast into an optimization problem for the functional 



J(E-T) = Tr{Op(T)}- 



X 



\E(t)f 



-dt - L 



(33) 



2 \J sin 2 (nt/T)' 

where O is the projection operator onto the target state, A is the penalty factor for 
strong fields. In this work A is treated as a Lagrangian multiplier such as to keep the 
integrated intensity close to Iq. To simplify matters we will not employ the HE approach 
at this point but resort to the Markovian approximation, which leads to a set of two 
coupled equations. 

Optimizing the functional J with respect to the field one gets the expression 



E(t) 



Tr{a(t)[/i,p(t)]}sin 2 (7rt/T), 



(34) 



where a(t) is an auxiliary operator propagating backward in time, starting from O at 
t = T with the following differential equation 



lh dt a 



-E(t)[p,a} + ihTZa, 



N 



rn 
N 



h m , a 



1,0- 



cr 



^(tot) 



a 



j(tot) 



(35) 



Eqs. (31) and (35) are nonlinear equations coupled through Eq. (34), which can be 



solved iteratively. 
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Figure 1. Level scheme for eigenstates of the consider model systems NA and 
NB. Black arrows mark the dominant transitions, grey arrows the next to dominant 
transitions starting from |li). 



The fields obtained by the OCT equations will be characterized by their XFROG 
(cross-correlation frequency-resolved optical gating [48] ) trace defined as 



j dt'E(t')G(t' - t)e~ iojt ' 



(36) 



where G(t) is the gate function having the form 

G(t) = J(erf[2(t + r/2)/A] + erf[2(* - r/2)/A]) . 



(37) 



In the above equation r is the width of the gate and A is the width of the shoulder. 



3. Results 



3.1. Parameters 

We will apply the formalism outlined in Section [2] to a model mimicking the 
situation in perylene bisimide aggregates. In particular we consider the case 
of A^,A^'-Di(iV-(2-aminoethyl)-benzamide)-l,6,7,12-tetra(4-tert-butylphenoxy)-3,4:9,10- 
perylenebiscarboximide (PBI) which is known to form J-aggregates [49]. The 
photophysical properties of PBI monomers have been studied in Ref . [JT] • As far as the 
present model is concerned the following properties will be used: The So-Si transition 
energy is E me = 2.13 eV. The monomer S -Si transition dipole moment is 3.34 ea . 
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This transition is coupled to an effective vibrational mode of frequency 1415 
cm -1 with a Huang- Rhys factor of 0.44. Further we used n — 1 in Eq. (19). While 
the effective mode captured the vibrational side band observed in the experiment, it 
could not account for the general line broadening, which amounts to 1110 cm -1 thus 
indicating substantial coupling to a continuous distribution of low-frequency bath 
modes. Hence, for the local high-frequency mode we take in the spectral density, Eq. 
(26) Uq = 1415 cm -1 and 7/h = 0.44, together with a phenomenological broadening 
width Ah = 100 cm -1 . For the low-frequency heat bath consisting of intermolecular and 
solvent vibrations, a frequency cut-off of A L = 100 cm -1 is used in Eq. (28 ) together with 
a coupling strength of rj^ = 2.0, chosen such as to yield considerable line broadening. 
The simulation is carried out at a temperature of 300 K, i.e. we can use Ah = Al = 1 
in Eq. (30). Convergent results have been obtained for a hierarchy order of seven. 

The situation is more complicated for the monomeric excited state absorption 
(ESA). On the experimental side, the ESA spectrum has been obtained from a pump- 
probe spectrum, recorded after equilibration in the Si state. Since the procedure requires 
knowledge about ground state bleach and stimulated emission spectra at a detail which 
cannot be obtained, the extracted ESA spectrum may contain spurious features. On 
the theoretical side calculating highly excited states for a molecule as large as PBI is a 
challenge, in particular because standard time-dependent density functional theory will 
fail to describe transitions of double excitation character. In Ref. [41J we applied the 
multireference variant of density functional theory to a reduced PBI system. Combining 
experimental and theoretical data it can be stated that (i) the excited state responsible 
for the Si-S n transition is at n > 20. In other words, the internal conversion down 
to the Si state proceeds via a dense manifold of electronic states, whose microscopic 
description is out of reach. This suggests using a description in terms of an effective 
internal conversion rate as outlined in Section [2] For the actual value we have used an 
inverse rate of 100 fs. Within the spectral density model, Eq. (24), this amounts 
to choosing Aic=1000 cm -1 and 771c = 0.4. (ii) Since reliable information on the 
transition dipole moment /i m j are not available we will use the value obtained within the 
harmonic oscillator picture, i.e. = \/2jjL mfi [15]. (iii) Concerning the anharmonicity 
A m = E m j — 2E n ^ e we first note that in the range where E m j w 2E n ^ e two transitions 
have been observed. Since the magnitude of A m decides about the mixing between local 
and nonlocal double excitation states [15] we will consider two cases. In case A the 
anharmonicity is large, A m = —0.26 eV, whereas in case B it is small, A m = —0.04 eV. 

In Ref. [H] a PBI monomer has been investigated only. The Coulomb coupling 
between adjacent monomers in the PBI aggregate has been estimated from the 
experimentally observed absorption line shifts upon aggregation. It has also been 
calculated in Ref. [5D]. Below we will use the calculated value of J mn = —515 cm -1 and 
only nearest neighbor coupling will be considered. The coupling between excited state 
transitions is not known and we will again use the harmonic oscillator approximation 

giving J^nn = VZJmn- 

Aggregates of three different sizes (N =2, 3, and 4 sites) will be studied for cases A 
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Figure 2. Population dynamics in the field- free case after initial population of a 
specific two-exciton eigenstate. Shown is the total population of the two-exciton 
manifold in comparison with the monomer case. 
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case B 
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Figure 3. Population dynamics of the dimer in the field-free case after initial 
population of the highest state of the two-exciton band. Shown are the populations of 
the one- and two-exciton eigenstates. 



and B (notation NA/NB). The one- and two-exciton eigenstates are shown in Fig. [T] 
Clearly, for cases A representing a larger anharmonicity the locally double excited states 
are energetically separated from the delocalized double excitation. For instance, in the 
dimer we have for case A |2 X > = 0.654 |11) +0.654 |22) +0.378 |12), |2 3 ) = —0.268 |11> — 
0.268 1 22) + 0.926 1 12) , i.e. the lowest/highest state has a dominant local/nonlocal 
character, whereas for case B we have |2i) = 0.537 1 11) + 0.537 1 22) + 0.650 1 12) and 
|2 3 ) = -0.460 |11) - 0.460 |22) + 0.760 |12), i.e. local and nonlocal double excitation are 
mixed. Similar results are found for N = ?> and 4. 

3.2. Population Dynamics in Field-Free Case 

In Fig. [2] we compare the internal conversion dynamics for the different models as 
a function of the initial preparation of the two-exciton manifold. Also shown is the 
decay of the monomer LDE state, which is always faster than that for the aggregate. 
Furthermore, it is found that the decay in case B is always slower as compared to case 
A, with the difference becoming more apparent upon increasing the aggregate size. This 
observation can be traced to the fact, that the mixing of LDE and NDE states in case B 
leads to a different intra-band dynamics and in particular to a slow-down of the inter- 
band dynamics as can be seen from Fig. [3j From this figure we notice that even though 
the decay of the initial state population is faster in case B, the population gets trapped 
for a longer time in state |2i), i.e. at the bottom of the two-exciton band. In case A the 
state |2x) is of more local character as compared with case B and therefore inter-band 
relaxation is faster. Similar arguments apply to the aggregates with N = ?> and 4. 

So far we have presented results from the HE simulation. In Fig. [4] we show a 
comparison between HE and Markovian dynamics. In case B one hardly notices any 
difference in the dynamics, whereas case A shows considerable difference between the 
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Figure 4. Population dynamics in field-free case after initial population of the highest 
and lowest state of the two-exciton band (top row: dimer, bottom row: tetramer). 
Shown are the total populations of the two-exciton bands for the HE and the Markovian 
case. 



HE and the Markovian dynamics if the system is prepared in the uppermost two-exciton 
state. In order to scrutinize this effect we have plotted the population dynamics of the 
two-exciton eigenstates in Fig. [5} First, we notice that in general the populations from 
the HE are showing an oscillatory behavior, which is not present in the Markovian 
approximation. Despite this difference in details, the total two-exciton population is 
rather similar in the two limits. Only for case A and for initial preparation in state 1 23) 
we notice a marked deviation for the dimer considered in Fig. [5j In the HE case state 
|2 3 ) relaxes rapidly and states |2 2 ) and |2i) become populated before the overall decay 
due to IC sets in. Notice that the energy gap between states |2 3 ) and |2 2 ) amounts to 
2518 cm -1 . The spectral density for the system-bath interaction is composed of a low- 
frequency part with cut-off at 100 cm -1 and a high-frequency contribution at 1415 cm -1 . 
Hence relaxation within second order perturbation theory as implied in the Markovian, 
i.e. truncated HE, approach will be very inefficient due to the smallness of the spectral 
density. The HE, on the other hand, accounts for higher order effects and yields a rapid 
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Figure 5. Population dynamics of two-exciton eigenstates for dimer. Shown are the 
population in each state starting from the lowest (top row) and highest (bottom row) 
state in the two-exciton band. 



energy relaxation. Notice that in case of model B this energy gap is only 1203 cm -1 , i.e. 
within the frequency range covered by the spectral density. According to the argument 
given one would not expect a marked difference between HE and Markovian dynamics 
in case B, what is in line with the results shown in Fig. |4| Inspecting Fig. [T]we note 
that a similar energy gap exists also for larger aggregates what leads to a corresponding 
behavior of the population dynamics in Fig. |4| However, upon further increase of the 
aggregate size this gap closes and it can be expected that HE and Markovian dynamics 
behave rather similar at least from the total population point of view. In passing we note 
that ultrafast inter-band relaxation suppresses the effect of excitonic-polaron formation, 
which, however, will play a role in the one-exciton manifold (see, Ref. [39]). 

3.3. Laser Control of Two-Exciton Dynamics 

In the following we will investigate the possibility to control two-exciton dynamics with 
ultrashort shaped laser pulses. In particular we will ask the question whether one can 
transiently suppress EEA. Since EEA is a local process one might argue that a two- 
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exciton wave packet, which is prepared in a way such that the contribution coming 
from LDE is minimized will lead to a slow down of EEA since the latter is mediated by 
intra-band relaxation processes mixing LDE and NDE states. For this purpose we will 
use the target state \1N) within the OCT scheme, i.e. the state corresponding to the 
situation where the two NDEs have the largest separation in real space. The fastness of 
the IC process puts some restriction on the pulse duration which has been set to T = 50 
fs. As mentioned in Sec. 2.5 the field-optimization will be performed in the Markovian 
limit, starting with a very broad band pulse. The field is then used to generate the 
dynamics using the HE approach. 

In Fig. [6] we show the optimized fields (XFROG) for cases A and B together with 
the population dynamics of target states as well as the exciton eigenstates for the dimer 
model. Since dipole transitions are possible between adjacent exciton bands only, the 
dynamics necessarily involves a two-step process. First, the one-exciton band is excited 
(transition frequency hcu\ lt0 = 2.06 eV) and this process is followed by a one- to two- 
exciton transition. In case A this way a superposition of states |2i) and 1 23) is prepared 
(transition frequencies hu>2 1 ,i 1 = 1-88 eV, ^2 3 ,ii = 2.06 eV), which, however, rapidly 
decays due to intra- and inter-band relaxation processes. For case B, where the LDE 
and NDE states are strongly mixed, the compromise found by the OCT equations has 
been to prepare just state |2i) which has a 42 % overlap with the target state. Although 
state |2 3 ) would have had a larger overlap, its transition dipole matrix element is about a 
factor of 13 smaller for that case (transition frequencies: hu 2lt i 1 = 2.04 eV, ^2 3 ,ii = 2.30 
eV). Similar to the scenario of the field free dynamics in Fig. |3j the overall decay of the 
two-exciton population is slower in case B as compared to case A. Finally, we notice that 
the maxima of the XFROG traces do not match with the bare exciton transitions in 
all cases, since superposition states are prepared by the broadband excitation. Further, 
exciton-vibrational coupling and the dynamic Stark shift will modify the bare excitonic 
resonances. 

One might argue that in the dimer EEA will always be very effective since the NDE 
are localized at neighboring sites. This situation might change for larger aggregates. As 
an example we consider the tetramer model in Fig. [7} The convergence of the OCT 
algorithm is rather slow in this case and for the given constraints only a small population 
of the target state can be reached. In both cases the pulse initially excites the lowest 
transition of the one-exciton band (ftco>i lj0 = 2.03 eV). For case A the subsequent one- 
to two-exciton transition populates mainly state |2s) (hw2 5 ,i 1 = 2.12 eV) which has 
the largest overlap with the target state (32 %). In case B state |2i) is dominantly 
populated {fou)2 1 ,i l = 2.01 eV) although it has only a small overlap with the target state 
(3%). However, as compared with state |2 4 ), which has a 37 % overlap with the target 
state, the transition dipole moment to state |2i) is larger by a factor of about 5. Overall, 
the comparison between cases A and B resembles again Fig. |6j 

Comparing the dimer and tetramer cases one notices that the initial hypothesis 
that the longer the aggregate the longer the time scale during which one can maintain 
a two-exciton population does not hold in general. Instead the difference between cases 




Figure 6. Laser-driven dynamics of the two dimer models. The target state has been 
1 12) and J in Eq. @ was set to 0.005 a.u.. Upper row: XFROG traces (Eq. @) 
of optimized field (20 (A) and 135 (B) iterations, parameters for XFROG are r = 5 
fs, A = 1 fs). Contours are the equi-interval plots with 5% increment of the peak 
height. The maximum field amplitude is 2.6 GV/m for case A and 3.1 GV/m for 
case B. The dashed line denotes the one-exciton resonance, the solid lines the one- to 
two-exciton resonances (colors as in Fig. []}. Middle row: Population dynamics of one- 
and two-exciton eigenstates. Lower row: Dynamics of target state and total one- and 
two-exciton population. 
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Figure 7. Laser-driven dynamics of tetramer models. Upper row: XFROG traces 
(Eq. (36)) of optimized field (38 (A) and 30 (B) iterations, parameters for XFROG 
are r = 5 fs, A = 1 fs). Contours are the equi- interval plots with 5% increment of the 
peak height. The maximum field amplitude is 2.3 GV/m for case A and 2.4 GV/m for 
case B. The dashed line denotes the one-exciton resonance, the solid lines the one- to 
two-exciton resonances (colors as in Fig. [I]). Lower row: Dynamics of target state and 
total one- and two-exciton population. 



A and B points to the importance of mixing between LDE and NDE states within the 
two-exciton band. However, focussing on case B only there is the anticipated increase 
of the time-scale of transient two-exciton state population. In this case the optimized 
pulse populates the lowest state of the two-exciton manifold whose overlap with the 
target state decreases with increasing system sizes due to the "dilution" of the zero- 
order excitation states. For the same reason, however, the overlap of state \2\) with 
LDE states decreases yielding a longer time scale for the two-exciton decay. 

Next we comment on the use of pulses derived by using the OCT equations and 
Markovian dynamics within the HE scheme. Inspection of the different cases shows 
that the resulting population of the target state is comparable in the two cases, i.e. this 
procedure does not lead to a degradation of control for the given pulse shape. Needless 
to say that in line with the discussion of Fig. [4] the relaxation dynamics after the pulse 
is over differs. Of course, one might argue that the combination of OCT equations and 
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Figure 8. Comparison of laser-driven dynamics for the dimer models using the pulses 
obtained from the OCT equations (red) and a fit of these pulses to two Gaussians 
(2G, black). The lower panels show the resulting target and two-exciton manifold 
populations. 



HE might result in different pulses. However, in view of the difficulties arising from 
the short life time of the two-exciton population this will have little relevance for the 
present discussion. 

We need to emphasize that our model is limited to the two-exciton space. In 
principle for the given field intensities it is likely that higher exciton manifolds could 
be excited as well (for a study of intensity-dependent multi-exciton dynamics within a 
Bloch model, see Ref. [52]). The expected rapid decay of these states would lead to 
an additional channel for populating the two-exciton manifold, what could in principle 
influence the details of the dynamics. Needless to say, that this effect could be suppressed 
by reducing the field intensity, but at the expense of the population of the two-exciton 
manifold. 

Finally, the question arises whether all the details of the control pulses play a role 
or in other words how robust is the achieved population control and can the pulse be 
simplified. Exemplarily, in Fig. [8] we show a comparison of the dynamics obtained for 
the OCT derived pulses (cf. Fig. [6]) and their fits to two Gaussian shaped pulses for the 
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dimer model. For case A we observe a strong sensitivity, i.e. the two-exciton population 
drops down considerably when using two Gaussian pulses. Apparently, the preparation 
of the |2i) and |2 3 ) superposition state depends on the details of the control field. In 
contrast, in case B where dominantly the state |2i) is prepared a simple Gaussian-shaped 
field is sufficient to achieve a control comparable to the optimized pulse. 



4. Summary 

Two-exciton dynamics has been studied on the basis of a non-perturbative and non- 
Markovian hierarchy equation approach. The dissipative dynamics has been combined 
with optimal control theory for obtaining laser fields designed such as to trigger long- 
lived two-exciton state populations. Specific results have been obtained for short 
aggregates made of J-aggregate forming perylene bisimide dyes for which Frenkel exciton 
parameters are available jH] . Establishing a two-exciton population one has to compete 
with the very efficient internal conversion, which is a consequence of the breakdown of 
the Born-Oppenheimer approximation. For the monomer this process takes place on a 
time scale of about 100 fs, which therefore sets the upper bound for laser control. The 
fact that two-exciton populations can be maintained on a longer time scale is due to the 
mixing between local and non-local double excitations of the aggregate; the latter do 
not decay on an ultrafast time scale. The considered dyes support two possible double 
excitation states in the relevant energy range. Therefore, we considered two scenarios 
corresponding to small and substantial mixing between the two type of aggregate 
excitations. It turned out that the case of strong mixing allows for maintaining a two- 
exciton population on a longer time scale (in the present cases about 1 ps as compared 
with 0.5 ps for the weak mixing case). This effect should be observable in a pump-probe 
experiment. At this point it should be noted that two-exciton populations in aggregates 
made of organic dyes like PBI have been observed for the first time only very recently 
[53] . what demonstrates feasibility of preparation and spectral identification. 

The present investigation highlights the importance of zero-order state mixing 
within the two-exciton band for the inter-band transitions. The more diluted the zero- 
order states the more the decay is slowed down. However, this situation might change 
if static disorder is taken into account which will lead to a localization of the two 
excitations on different parts of the aggregate. 



Appendix: Derivation of the Hierarchy Equations of Motion 



For the Caldeira-Leggett model of dissipation, Eq. (19), with Hs-b = f(s)g(b), one 
can employ a stochastic procedure to derive the equation of the motion for the reduced 
density matrix of the relevant system p if the whole system is prepared as a factorized 
state, i.e. ptot(0) = p(0) exp(— (3H B ) /Tr[exp(— 0Hb)] J3TJ, [38] . Within this approach, the 
influence of the bath is completely characterized by its response function and its effect 
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is described by a random force g(t) in the Ito stochastic differential equation 

ihdp = [H s + g(t)f(s),p s ]dt + Vh/2 [f(s),p] dw 1>t + iVh/2 {f(s),p} dw* 2>t , (A.l) 

where w\ tt , w 2jt are two independent complex-valued Wiener processes and w* t , w 2t are 
their complex conjugates. In Eq. (A.l) g(t), which plays the same role as the influence 
functional in the path integral treatment, can be regarded as a stochastic field fully 
characterizing the influence of the environment. It is defined as 



(A.2) 



g(t) = Vh J ^a R (t — t')dw\ t , + oti{t — t')dw 2 ^. 
For simplicity, here it is assumed that the response function of the bath is 

2 

a(f) = X>*<T n **, (A.3) 



k=l 



where Vt\ = Q 2 ( c ^- Section 2.3.2). Note that the extension to multiple exponentials is 
straightforward. For a more general form of the response function consult Ref. 
For the case of Eq. (A.3) the stochastic force can be written as 

-dw 2 ,r) 



g(t) = Vh / exp[— Qi(t — T)](ai ]+ dw\ T + a±- 
Jo 

+Vh / exp[-Q 2 (t - r)}{a 2 - + dw* 1 + a 2 ._dw 2 ) 
Jo 



= ^(t)+^ 2 (t). (A.4) 

Here a k -± are defined as oti ;+ = (b x + b 2 )/2, a 1; _ = -i(bi - b 2 )/2, a 2 - + = (b* + b 2 )/2, 
and a 2 -- = -i{-b\ + b 2 )/2. 

Next we introduce the auxiliary reduced density matrices 

p m ,n(t) = J^ESELm mmum)} , (a.5) 

ym!n!(c[ H) ) m+ " 

where M{ - ■ •} denotes the ensemble average with respect to the noise and c H ^ = \f\b\b 2 \. 
Notice that using a scaling like in the above equation has been suggested by YiJing Yan 
and co-workers [36, 51J. The present prefactor is slightly different from that suggested 
in Refs. [361 EI]; where the numerator was set equal to one. This choice will keep 
the terms in the same order size consistent in cases where the decay constants in the 
response functions are rather different. 

Carrying out the stochastic average for the random variables in Eq. (A.l), 
elementary stochastic calculus yields the differential equation for p mn (t) 

ihdp m>n (t) / dt = -ih^rnQx + nft 2 )p m>n (t) + [H 8 , p m>n (t)] 
+^ 1 )- 1 V / (m + l)4 H) [/(s),p m+1 ,„,(t)] 

+^ 2 )- i v / (n + i)4 H) [/(s,p m , n+1 (t)] 

+Mt(n 1 )^/m/c n) (a 1;+ [f(s) ,Pm-l,n{t)} — Q!i ; _{/(s), pm-l, n {t)} 



+hft(n 2 )^n/ C W(a 2 .,4f(s),p m , n - 1 (t)] - a 2 ^{f(s),p m>n ^(t)}) (A.6) 
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This equation needs to be solved for the initial conditions poo(0) = p(0) and p m ,„(0) = 
(m ^ orn / 0). Application to the present system-bath model yields the equations 
of motion, Eq. (30 ). 
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